
set more off
set matsize 11000

cd \Data

********************************************************************************
** Create continuous temperatures values to search for spline minimums
// Minimum avg. daily temperature is 31.71 (this will serve as the lower bound)
// Maximum avg. daily temperature is 88.88 (this will serve as the upper bound)

set obs 5718
gen temp_avg = _n
replace temp_avg = 31.70 + temp_avg/100
mkspline temp = temp_avg, cubic knots(52 62 72)
save Temp/min_temp_search.dta, replace


********************************************************************************
** Estimate premise-level restricted cubic spline models (knots at 52, 62, 72)
** Note: the output for each premise will be a row-vector with 5 elements:
** hh_`id' = [premise_id; predicted aggregate cooling kwh during 2012 and 2013; actual consumption during 2012 and 2013; spline minimum height; temperature at which spline minimum is reached]
********************************************************************************

** Note: the cons_data.dta dataset is split up into 16 subsets

forvalues j = 1/15 {

	local id_min = 10000 * (`j' - 1) + 1
	local id_max = 10000 * `j'
	
	use cons_data.dta, clear
	keep if premise_id >= `id_min' & premise_id <= `id_max'
	mkspline temp = temp_avg, cubic knots(52 62 72)
	save Temp\cons_data_`j'.dta, replace

}

use cons_data.dta, clear
keep if premise_id >= 150001 & premise_id <= 158112
mkspline temp = temp_avg, cubic knots(52 62 72)
save Temp\cons_data_16.dta, replace
	

forvalues j = 1/15 {

	local id_min = 10000 * (`j' - 1) + 1
	local id_max = 10000 * `j'
	
	forvalues i = `id_min'/`id_max' {
	
		use Temp/cons_data_`j'.dta, clear
		keep if premise_id == `i'
		quietly: reg kwh_daily temp1 temp2
		append using Temp/min_temp_search.dta
		predict fitted
		quietly: sum fitted if kwh_daily == .
		scalar define min_kwh = r(min)
		quietly: sum temp_avg if fitted == min_kwh & kwh_daily == .
		scalar define min_temp = r(mean)
	
		drop if kwh_daily == .
		gen cooling_kwh = 0
		replace cooling_kwh = fitted - min_kwh if temp_avg > min_temp
		collapse (sum) cooling_kwh kwh_daily, by(premise_id)
	
		mkmat premise_id cooling_kwh kwh_daily
		matrix define hh_`i' = (premise_id[1,1], cooling_kwh[1,1], kwh_daily, min_kwh, min_temp)
	
	}
	
	matrix define subset_`j' = (hh_`id_min')
	local k = 10000 * (`j' - 1) + 2
	forvalues i = `k'/`id_max' {
		matrix define subset_`j' = (subset_`j' \ hh_`i')
	}
}


forvalues i = 150001/158112 {
	
	use Temp/cons_data_16.dta, clear
	keep if premise_id == `i'
	quietly: reg kwh_daily temp1 temp2
	append using Temp/min_temp_search.dta
	predict fitted
	quietly: sum fitted if kwh_daily == .
	scalar define min_kwh = r(min)
	quietly: sum temp_avg if fitted == min_kwh & kwh_daily == .
	scalar define min_temp = r(mean)
	
	drop if kwh_daily == .
	gen cooling_kwh = 0
	replace cooling_kwh = fitted - min_kwh if temp_avg > min_temp
	collapse (sum) cooling_kwh kwh_daily, by(premise_id)
	
	mkmat premise_id cooling_kwh kwh_daily
	matrix define hh_`i' = (premise_id[1,1], cooling_kwh[1,1], kwh_daily, min_kwh, min_temp)
	
	}
	
matrix define subset_16 = (hh_150001)
forvalues i = 150002/158112 {
	matrix define subset_16 = (subset_16 \ hh_`i')
}



clear
forvalues j = 1/16 {

	svmat subset_`j'
	rename subset_`j'1 premise_id
	rename subset_`j'2 cooling_kwh
	rename subset_`j'3 total_kwh
	rename subset_`j'4 min_kwh
	rename subset_`j'4 min_temp
	save Temp/hh_results_`j'.dta, replace
	clear
	
}
	
use Temp/hh_results_1.dta, clear
append using Temp/hh_results_2.dta
append using Temp/hh_results_3.dta
append using Temp/hh_results_4.dta
append using Temp/hh_results_5.dta
append using Temp/hh_results_6.dta
append using Temp/hh_results_7.dta
append using Temp/hh_results_8.dta
append using Temp/hh_results_9.dta
append using Temp/hh_results_10.dta
append using Temp/hh_results_11.dta
append using Temp/hh_results_12.dta
append using Temp/hh_results_13.dta
append using Temp/hh_results_14.dta
append using Temp/hh_results_15.dta
append using Temp/hh_results_16.dta
save Temp/HH_Cooling_Estimates.dta, replace



